%% master_counterfactuals.m


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Produced for JdG,AL,IM by Christopher Evans at UPF 
%
% Date: February 2019
%
% Purpose: This file runs the counterfactual exercises after the
% main_baseline.m has been run. Change the shocks to experiment with the
% counterfactuals.
% 
% Based on one file: Step_4_Algorithm_counterfactual.m, which uses new
% baseline input data, MAT/Step_4_setup_`alphaT''rhoT'_`shockT'.mat
% as initial level variables, and then shocks the system and computes hat
% variables for wages, prices, exports, etc.
%
% NEW: for CES setup. Have streamline file that has to be called up and
%                     various variables that have to be used.
% Output: MAT/Step_4_results_`alphaT''rhoT'`lambdaT'`etaT'_`shockT'.mat,
%         MAT/Step_5_results_`alphaT''rhoT'`lambdaT'`etaT'_`shockT'.mat
%         MAT/Elasticity_output_`alphaT''rhoT'`lambdaT'`etaT'_`shockT'.mat
%
% Last Updated: 07/09/2019 by JDG/IM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global alphaT rhoT lambdaT etaT shockT rho sigma sigmaT lambda eta ...
    tol maxit Pi_hat_n D_hat_n Xi_hat_mnjf Xi_hat_mnj a_hat_f a_hat ...
    varphi sL_n sPi_n sD_n N J FI_J france oneminus_alpha_vector ...
    labour_share_PWT countrysecD firmsecD_sorted start start_sorted ...
    sum_by_country_dummy ISO alpha_WIOD_j alpha_WIOD_francej sector_algorithm 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% alphaT: 
%
% 1 = PWT labor share for all countries and Labor share for French firms 
%     are at the COUNTRY level.
%
% 2 = PWT labor share for all countries and WIOD scaled labor share for 
%     French firms. Labor share for French firms are at the FIRM level.
%
% 3 = WIOD labor share for all SECTORS with the labor share at the SECTOR 
%     level for French firms.
%
% 4 = WIOD labor share for all SECTORS with the labor share at the FIRM 
%     level for French firms 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%alphaT = 4;

%% Loading in the previous data to get some useful variables

fname = strcat('MAT/Step_2_firm_data_alpha',num2str(alphaT),'.mat');
load(fname,'FI_J', 'firmsecD_sorted', 'start_sorted', 'alphaT');
load('MAT/WIOD.mat','countrysecD','france','J','N','start','sum_by_country_dummy');

%% SET PARAMETERS VALUES

% Frisch elasticity for GHH preferences

%varphi=3;

% share of labour in total final consumption (need to make assumption or collect data)

sL_n = ones(1,N)*0.8; 

% share of profits in total final consumption (should have this data)

sPi_n= ones(1,N)*0.15; 

% share of trade deficit in total final consumption (should have this data)

sD_n= ones(1,N)*0.05; 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% rhoT:
%
% Value for rho. Will apply to J number of sectors to create a vector of
% homogeneous rhos across sectors.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%rhoT = '3';

%rho = 3*ones(J,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Value for sigma
%
% This is the elasticity of substitution across sector consumption. Given
% we do not adjust this anymore after running C-D production we set to the
% value we think is most reasonable: 1.5. There is therefore no longer a
% "type" variable.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%sigma = 1.5*ones(J,1);
%sigmaT= '1.5';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% lambdaT:
%
% Value for lambda. We might experiment on this so have a type.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%lambdaT = '1.000001';

%lambda = 1.000001;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% etaT:
%
% Value for eta. We might experiment on this so have a type.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%etaT = '1.000001';

%eta = 1.000001;

%% Baseline: no shocks

% Pi_hat: Change in profits per country
    
Pi_hat_n=ones(1,N);
    
% D_hat: Change in trade deficit
    
D_hat_n = ones(1,N);
    
% Taste shock - for France only
    
Xi_hat_mnjf=ones(FI_J,N);
    
%Taste shock - for sectors only (countries except France)
    
Xi_hat_mnj= ones(N*J,N);
    
% Productivity shock, a is firm specific - France only
    
a_hat_f= ones(FI_J,1);
    
% Productivity shock, a is for sectors - countries except France
    
a_hat = ones(J*N,1);

for qq = [2,4]
    %1:4
    
%% Feeding in SHOCKS
    
    if qq == 1
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Shock 1:
    %
    % A +10% productivity shock to all firms for Germany. Germany is the 10th
    % country in the WIOD.
    % UNIQUE IDENTIFIER (shockT) keeps track of counterfactual exercise
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    shockT = 'DEUprod_p10'
    
    a_hat = ones(J*N,1);
    a_hat(start(10):start(10+1)-1)=1/1.1;
    Xi_hat_mnjf=ones(FI_J,N);
    
    elseif qq == 2 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Shock 2:
    %
    % A +10% productivity shock to all firms to all countries except France.
    % UNIQUE IDENTIFIER (shockT) keeps track of counterfactual exercise
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    shockT = 'WORLDprod_p10'
    
    a_hat = ones(J*N,1);
    a_hat=a_hat/1.1;
    a_hat(start(15):start(15+1)-1)=1;   % Resets France to no shock
    Xi_hat_mnjf=ones(FI_J,N);
    
    elseif qq == 3
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Shock 3:
    %
    % A +10% preference demand to all firms from Germany. Germany is the 10th
    % country in the WIOD.
    % UNIQUE IDENTIFIER (shockT) keeps track of counterfactual exercise
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    shockT = 'DEUpref_p10'
    
     Xi_hat_mnjf=ones(FI_J,N);
    Xi_hat_mnjf(:,10) = 1.1^(str2num(rhoT)-1);
    a_hat = ones(J*N,1);
    
    else
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Shock 4:
    %
    % A +10% preference demand to all firms from the World. France is the 15th
    % country in the WIOD.
    % UNIQUE IDENTIFIER (shockT) keeps track of counterfactual exercise
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    shockT = 'WORLDpref_p10'
    
    Xi_hat_mnjf = ones(FI_J,N)*1.1^(str2num(rhoT)-1);
    Xi_hat_mnjf(:,15) = 1;
    a_hat = ones(J*N,1);
    end
  
    %% Running the Algorithm
    
   disp('BEGINNING COUNTERFACTUAL EXERCISE')
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
    % Step_4_Algorithm: same algorithm as Step_4_Algorithm_baseline.m, but now
    % we are able to feed in the model the calculated t+1 varialbes as t,
    % starting at the new equilibrium with no wedges
        eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_4_setup_alpha',num2str(alphaT),...
    '_rho',rhoT,'_approx.mat')]); 
      Step_4_Algorithm_counterfactual_CES_approx
%     
%     % Save the counterfactual output so it can be used for post analysis
% 
      outputtable = [' X_mnj X_hat_mnj pi_mnjf1 pi_mnj1 ' ...
         'pi_c_mnj1 pi_c_nj1 pi_M1 pi_l1 pi_M_f1 pi_l_f1 '...
          'PC_hat_n1 PC_n P_hat_n P_hat_nj P_hat_mnj w_hat ' ...
          'a_hat a_hat_f Xi_hat_mnj Xi_hat_mnjf '...
          'final_goods1 intermediates1'] ;
% %  
      eval([strcat('save MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_4_results_alpha',num2str(alphaT),...
         '_rho',rhoT,'_',shockT,'_approx.mat',outputtable,' -v7.3;')]);   
%      
   display('STEP 4 COUNTERFACTUAL COMPLETED SUCCESFULLY')
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    % Step_5_counterfactual_output: Allows the user to easily analyse the
    %      counterfactual exercises produced before. Runs the CounterOutput.m
    %      function after specifying alpha type (labor share), gamma
    %      calibration and the shock scheme to be analysed.

                %% Load data

            
            % Time 0 data

            eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,...
    '_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_4_setup_alpha',num2str(alphaT),...
     '_rho',rhoT,'_approx.mat')]);

            % Time 1 data

            eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,...
     '_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_4_results_alpha',num2str(alphaT),...
     '_rho',rhoT,'_',shockT,'_approx.mat')]);
 
     Step_5_counterfactual_output_CES_approx
    
    % Save the counterfactual output so it can be used for post analysis

      outputtable = [' VA_hat_n VA_hat_mj VA_hat_mnj VA_hat_fj VA_hat_fnj '...
          'RGDP_def_hat_n RGDP_def_hat_fj RGDP_cpi_hat_n RGDP_cpi_hat_fj '...
          'RGDP_doubledef_hat_n RGDP_doubledef_hat_fj '...
          'real_wage_cpi_hat real_wage_def_hat real_wage_doubledef_hat '...
          'real_wage_income_cpi_hat_n real_wage_income_def_hat_n real_wage_income_doubledef_hat_n '...
          'TFP_hat_def_n TFP_hat_cpi_n TFP_hat_doubledef_n TFP_hat_cpi_france_j '...
          'imports_over_GDP_all imports_over_GDP_france_sector '...
          'exports_over_GDP_all exports_over_GDP_france_sector '...
          'GDP_deflator_n DoubleDeflation_deflator_n P_hat_n VA_fj_0 VA_mj_0'] ;
      %'TFP_hat_def_france_j TFP_hat_doubledef_france_j '...
%  
     eval([strcat('save MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_5_output_alpha',num2str(alphaT),...
         '_rho',rhoT,'_',shockT,'_approx.mat',outputtable,' -v7.3;')]); 
  
%     
     display('STEP 5 COMPLETED SUCCESFULLY')
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Step_6_Elastisticity_output: perform covariance decomposition at firm and
    %        sector level

    % Load files: only needed if run separately from Step 5

    eval([strcat('load MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_5_output_alpha',num2str(alphaT),...
         '_rho',rhoT,'_',shockT,'_approx.mat;')]); 

    % Shock size
    
   delta = 0.1;
    
   Step_6_Elasticity_output_CES

       outputtable = [' e_Y_FRf_cpi m_e_f_cpi cov_e_f_cpi check_f_cpi '...
           'e_Y_FRs_cpi m_e_mjFRs_cpi cov_e_mjFRs_cpi check_mjFRs_cpi '...
           'e_Y_s_cpi m_e_mj_cpi cov_e_mj_cpi check_mj_cpi '...
           'e_Y_FRf_def m_e_f_def cov_e_f_def check_f_def '...
           'e_Y_FRs_def m_e_mjFRs_def cov_e_mjFRs_def check_mjFRs_def '...
           'e_Y_s_def m_e_mj_def cov_e_mj_def check_mj_def '...
           'e_Y_FRf_doubledef m_e_f_doubledef cov_e_f_doubledef check_f_doubledef '...
           'e_Y_FRs_doubledef m_e_mjFRs_doubledef cov_e_mjFRs_doubledef check_mjFRs_doubledef '...
           'e_Y_s_doubledef m_e_mj_doubledef cov_e_mj_doubledef check_mj_doubledef'] ;
%  
     eval([strcat('save MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_6_output_alpha',num2str(alphaT),...
         '_rho',rhoT,'_',shockT,'_approx.mat',outputtable,' -v7.3;')]); 

     display('STEP 6 COMPLETED SUCCESFULLY')
    
    disp('COUNTERFACTUAL EXERCISE COMPLETED SUCCESFULLY')
end